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I. Introduction 


The main objective of this thesis is to identify dynamical models of the Hydraulic 
Manipulator Test Bed (HMTB). In particular, system identification techniques will be 
used to identify the joint dynamics and to validate the correctness of the HMTB models. 
Though dynamic model verification has been studied and performed for the DOSS flight 
manipulator, dynamic system identification for the hydraulic kinematically-equivalent 
ground-based DOSS manipulator located in the hydraulic manipulator test bed (HMTB) 
facility at the NASA Langley Research Center has not been studied in detail. This thesis 
will describe, apply, and compare system identification techniques for three joints 
(shoulder yaw, shoulder pitch, and elbow pitch) of the seven DOF hydraulic manipulator 
for the purpose of obtaining an adequate dynamic model of HMTB during insertion of the 
remote power controller module ORU. 

To perform the identification, a series of single-input, single-output (SISO) and 
multi-input, multi-output (MIMO) experiments will be performed. Nonparametric and 
parametric identification techniques will be explored in order to develop representative 
models of the selected joints. The identified SISO model estimates will be validated. The 
best performing models will be used for a decoupled multivariable state-space model. It 
should be noted that each identified model represents an open-loop representation of the 
closed-loop implementation for each joint. It is not the purpose of this thesis to determine 
the effective inertia or the effective damping coefficients for the HMTB links. The 
manipulator is localized about a representative space station orbital replacement unit 
(ORU) exchange task allowing the use of linear system identification methods. The 
parametric models will be compared to determine the best dynamic model for performing 
the ORU task. 

System identification techniques have been applied in many different fields. The 
purpose of the identified models in this thesis is to use them in a control application. The 
thesis concludes by proposing a model reference control system to aid in astronaut ground 
tests. This approach would allow the identified models to mimic on-orbit dynamic 
characteristics of the actual flight manipulator thus providing astronauts with realistic on- 
orbit responses to perform space station tasks in a ground-based environment. 

The process of system identification starts by performing an identification 
experiment, that is, exciting the system using some sort of input signal and observing the 
output over a time interval [9], Once the experimental data is recorded, parametric or 
nonparametric analysis can be performed. In nonparametric analysis, a system’s transfer 
function, impulse response, or step response is extracted from the experimental data in 
order to determine transient or frequency response characteristics of the system. This 
method, however, is often sensitive to noise and usually does not give very accurate 
results [9], In parametric analysis, the recorded input and output sequences are fitted to a 
parametric model. This process begins by determining an appropriate model form. Next, 
some statistically based method is used to estimate the unknown parameters of the model. 
The model is then tested or validated to determine if it appropriately represents the 
dynamic system. 

The remainder of this chapter provides historical background of the DOSS 
manipulator, the Hydraulic Manipulator Test Bed (HMTB) housed at the NASA Langley 
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Research Center, and the orbital replacement unit hardware used by the manipulator. 
Most of this information has not been published before. The chapter concludes by 
providing a literature search on system identification techniques used in this thesis. 

Chapter II will describe the overall experiment design process developed 
specifically for the hydraulic manipulator test bed (HMTB). As a precursor to parametric 
identification. Chapter III will describe the application of nonparametric methods used to 
extract characteristics of the unknown joints. Parametric model estimation techniques 
primarily used for control system identification will be applied in Chapter IV. In this 
technique, transfer function models describing each joint and its associated disturbances 
are analyzed to yield an adequate state-space model approximation. The second 
parametric technique, used primarily in modal system identification, will be employed in 
Chapter V. This technique uses a minimum realization algorithm to determine a model 
with the smallest state-space dimension among all realizable systems. Comparisons of the 
parametric models will be shown in Chapter VI. Chapter VII concludes the thesis by 
providing suggestions for future work. A model reference control system is proposed to 
provide astronauts with realistic on-orbit responses to perform space station tasks on the 
ground. 

Matlab menu-driven system identification software programs were developed for 
this project. One of the programs, a menu-driven script written for nonparametric and 
parametric evaluation of the input/output data using functions from the MA TLAB System 
Identification Toolbox. Another menu-driven program was used to identify models using 
the Observer/Kalman Filter Identification (OKID) technique, provided in the 
System/Observer/Controller Identification Toolbox (SOCIT). This last program script 
used several toolboxes to perform MIMO comparisons for identified models. 

1 Dexterous Orbital Servicing System (DOSS') Background 

In 1984 President Reagan directed the National Aeronautics and Space 
Administration (NASA) to build a space station. He invited allies of the United States to 
join in the challenge of creating a machine that could be manned and operated beyond the 
year 2000 [1]. Space Station Freedom shown in Figure 1 was the first major co-operative 
program of the governments of the U.S., Japan, the 10 nations of the European Space 
Agency (ESA), and Canada for the utilization and operation of a microgravity laboratory 
environment in space. Each government was responsible for furnishing specific user 
elements of Space Station Freedom. The United States through the direction of the 
National Aeronautics and Space Administration (NASA) was responsible for the design, 
development, and construction of the truss assembly infrastructure, the crew living 
quarters (US Habitat Module), and the US Laboratory Module. Japan would develop and 
assemble the Japanese Experiment Module (JEM). The European Space Agency (ESA) 
and its member states would develop their own Free-Flying Laboratory named Columbus 
and a polar platform. Canada’s responsibility involved providing the Mobile Servicing 
System (MSS), a complex robotic machine used to assemble, service, and maintain most 
of the station. The MSS’s major robotic components are the Space Station Remote 
Manipulator System (SSRMS) and the Special Purpose Dexterous Manipulator (SPDM) 
shown in Figure 2. 
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Figure 2. Canada’s SSRMS and SPDM working on Freedom’s truss. 


The U S. Congress also appropriated a portion of space station money for U.S. supported 
space station robotics [17], With these funds, NASA started development of the Flight 
Telerobotic Servicer (FTS), a dexterous manipulator shown in Figure 3, for use on both 
the Space Transportation System (STS) and the space station. After determining the 
requirements for the space servicing manipulator, NASA awarded Martin Marietta 
Astronautics Group (MMAG) a contract to design, construct, and test a flight deliverable 
FTS system. 
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Figure 3. NASA’s Flight Telerobotic Servicer (FTS). 

When the Advisory Committee on the Future of the U.S. Space Program convened 
by Vice President Dan Quayle issued a report in December 1990, NASA’s goals, 
programs and practices were altered [2], The NASA program that suffered the most 
devastating blow was Space Station Freedom. The report recommended that Space 
Station Freedom be utterly recast, reduced in both scale and complexity, which was a 
decree previously urged by Congress [3], Freedom’s new primary mission was in the life 
sciences, specifically the psychological and physiological effects of microgravity on 
humans. In the space station redesign process, many services and capabilities were 
reduced while others were halted indefinitely such as NASA’s Flight Telerobotic Servicer 
(FTS) project. With over 270 million dollars already invested in the development and 
fabrication of FTS robotic technologies, NASA initiated a project apart from the space 
station to capture the newly developed FTS technologies. 

When Canadian politicians pushed to withdraw from the leaner, redesigned space 
station, NASA and the Canadian Space Agency (CSA) began talks in March 1994 to 
develop a plan which would reduce Canadian space station costs and bolster space science 
cooperation [4], In the new plan, Canada would defer the Special Purpose Dexterous 
Manipulator (SPDM) which is a significant portion of Canada’s robotic contribution to the 
space station. With SPDM, the space station’s primary robotic resource, deferred, NASA 
decided to continue a stunted version of its robotic program using existing FTS 
technology to provide a robotic presence on the shuttle and the space station in the 
interim. 

The new robotic thrust called the Dexterous Orbiter Servicing System was based 
on the previous FTS designed by Martin Marietta. Initially, the new system was to 
provide robotic capabilities to space shuttle astronauts. Mission specialists would utilize 
and test the system in the shuttle cargo bay as a precursor to space station related tasks 
and procedures. The project was later termed the Dexterous Orbital Servicing System 
(DOSS) to service both the space shuttle and the redesigned space station. 

A ground based trainer system composed of a shuttle aft flight deck (AFD) 
mockup and one seven-degree-of-freedom hydraulic manipulator mounted on a multi- 
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purpose experiment support structure (MPESS) was designed as a form, fit, and 
functional laboratory version of the flight system [5], The kinematically equivalent 
hydraulic manipulator was developed by Western Space and Marine (WSM). The ground 
based trainer system is referred to as the Hydraulic Manipulator Test Bed (HMTB). 

2. Hydraulic Manipulator Test Bed fHMTB) 

The DOSS trainer built by Western Space and Marine was first located at Martin 
Marietta and then transferred to NASA’s Langley Research Center (LaRC). The trainer 
was placed in the Hydraulic Manipulator Test Bed (HMTB) facility shown in Figure 4. 
The HMTB facility includes a ground test dexterous manipulator driven by Ada flight 
prototype software and a shuttle aft flight deck (AFD) mockup. 



Figure 4. LaRC’s Hydraulic Manipulator Test Bed (HMTB). 

Specifications were developed to train the flight crew to operate the trainer system 
in accomplishing mission tasks, to operate in a 1-G environment, and to develop mission 
operation timelines. Layout of the trainer system in the HMTB facility, shown in Figure 5, 
was configured to provide the flight crew the same geometry, camera views, and lighting 
conditions that would exist during the actual flight for completion of space related tasks. 
Mission tasks include the installation and removal of space station truss members, the 
exchange of space station orbital replacement units (ORU), mating thermal utility 
connectors, and performing inspection tasks. 

The trainer manipulator consists of seven hydraulic rack and pinion actuators and 
their controlling valves integrated with structure to provide a seven degree of freedom 
(DOF) hydraulic manipulator [6] as shown in Figure 6. The hydraulic manipulator 
provides the same kinematics as the flight manipulator, that is, six controllable degrees of 
freedom (shoulder yaw and pitch, wrist pitch, yaw and roll, and elbow pitch) with a single 
indexed roll DOF at the shoulder. 
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Figure 6. The HMTB Hydraulic Manipulator. 


6 


, -v;. 


The aft flight deck (AFD) mockup, shown in Figure 8, is a replica of the actual 
AFD (depicted in Figure 7) and provides some of its functions. It provides the crew with 
an interface to control telerobot tasks and operations with or without a direct view of the 
worksite. Additional devices within the AFD mockup include a 2x3 DOF hand controller 
to teleoperate the trainer arm, closed circuit television monitors to display views of the 
hydraulic trainer, and an emergency shutdown (ESD) switch to turn off power to each arm 
servo while maintaining power to the trainer’s control computer. Figure 8 displays an 
internal view of the AFD mockup. 



Figure 7. The Aft Flight Deck (AFD). 



Figure 8. Internal View AFD mockup. 

The trainer control station is comprised of a trainer control station computer, 
power control panels, a 1553B bus interface, and various other communication devices. 
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The trainer control station computer provides software for control of training, and setting 
up training conditions. The 1553B bus monitors commands and responses between the 
control computer, the joint controllers, and the 2x3 DOF hand controller. 

According to FTS trainer specifications, the HMTB facility at NASA LaRC is an 
adequate, ground based version of the system to be used by shuttle flight crew members in 
accomplishing mission tasks. HMTB configuration provides the crew with the same 
geometry, similar camera views, identical manipulator kinematic configuration, the same 
software control, and lighting conditions that would exist during an actual space shuttle 
flight. 

3 , Remote Power Controller Module 

One of the primary mission tasks on the space station will be the maintenance of 
orbital replacement units (ORUs). With approximately 70 remote power controller 
module (RPCM) ORUs located on various port and starboard clusters of the space 
station, extravehicular activities (EVA) performed by the space station crew members 
would be difficult, impractical, and potentially hazardous [7], Robotic servicing of the 
RPCM by the DOSS system would minimize the EVA crew time and significantly increase 
crew safety. 

There are six types of RPCMs each varying in power capacity while maintaining 
identical physical dimensions. Each RPCM, (mockup shown in Figure 9), is responsible 
for regulating and distributing secondary power to critical space station components. 
Therefore, replacement of failed RPCMs by a dexterous manipulator would provide a 
crucial space station maintenance service. Successful completion of an RPCM exchange 
procedure lies in the ability of the dexterous manipulator to extract the failed ORU, to 
exchange the failed ORU with a new one, and to carefully insert the new ORU. For the 
purpose of this thesis, data from the RPCM ORU exchange task will be used to validate 
different dynamical models for the HMTB joints. 



Figure 9. Remote Power Controller Module (RPCM) ORU Mockup. 
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4. System Identification Techniques 

With a myriad of interrelated approaches, perspectives, methods, techniques and 
specializations, the field of system identification has widespread application in many areas 
such as communications, geophysical engineering, fault detection, pattern recognition, 
adaptive filtering, linear prediction, electric circuits and robotics. With this in mind, the 
literature search has been limited to system identification techniques for control purposes. 

For thirty years system identification has been an important discipline within the 
area of controls. With modem control methods requiring specific accuracies for 
mathematical models, the system identification benefit of improving an analytical model is 
of noted significance. Two complementary aspects of system identification are frequency 
domain and time domain identification. 

Frequency domain identification historically dominated system identification 
practice in control engineering prior to the 1960s [11]. Frequency domain identification 
emphasizing nonparametric identification methods has been used for stability, design and 
control purposes. In this thesis, transfer function, correlation, and spectral analysis 
techniques will be used for nonparametric identification of the HMTB joints. 

Time domain approaches emphasize parametric identification techniques for the 
system identification problem. The last two decades have seen a tremendous increase in 
the use of parametric time domain identification methods. This has been partly due to 
stricter accuracy requirements for mathematical control models as well as the increased 
availability of digital computers that can estimate system characteristics much faster than 
conventional frequency domain methods. The first parametric identification technique 
employed in this thesis will use various black-box transfer function model structures to 
determine parametric model estimates for the HMTB joints [10], The transfer function 
models, found within the MATLAB System Identification Toolbox by Lennart Ljung [10], 
use a prediction error method (PEM) to determine parameters for each black-box model. 
The PEM is a modification of the least squares (LS) method. 

The field of structures has used parameter identification techniques based on 
system realization theory. One such technique, Observer/Kalman Filter identification 
(OKID), will also be used to identify parametric models for the HMTB joints. This 
minimum realization approach to time domain system identification yields a model with the 
smallest state space dimension among a set of models having the same input-output 
relationship. Ho and Kalman [12] both developed minimum realization theory using 
Markov parameters which are simply pulse response functions. 

In this thesis, both frequency and time domain techniques will be used to extract 
and identify dynamic characteristics of the HMTB manipulator. The particular techniques 
used in this thesis will be discussed in more detail after the following chapter which will 
describe the experiment design process and setup. 
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n. Experiment Design 


This chapter will describe the overall experiment design process developed 
specifically for the hydraulic manipulator test bed (HMTB) at the NASA Langley 
Research Center. The experiment design has been modularly configured and developed 
within physical hardware limitations and temporal constraints. 

1 . Overall Design Process 

The steps of identifying dynamic models of the manipulator joints involve 
designing an experiment, selecting a model structure, choosing a criterion to fit, and 
devising a procedure to validate the chosen model. With the goal of obtaining a ‘good 
and reliable model estimate, Ljung [10] emphasizes the importance of the experiment 
design and the selection of its associated variables. Since a good model is not likely to be 
obtained from bad experiments, identification experiments should be selected to effectively 
characterize all the important modes of the system. This involves selecting persistently 
exciting (pe) input signals, that is, input signals which have strictly positive spectral 
density functions for all frequencies in the frequency band which is of interest for the 
intended application of the model. 

Figure 10 displays a pictorial representation of the experiment design segments 
starting from the generation of excitation signals to the extraction of useful joint space 
data for system identification and parameter estimation. Several physical devices and 
software applications were used in the experiment design process. Generation of the 
various waveform signals was performed on an IBM compatible 486 personal computer 
(PC). 



Waveform Excitation HMTB Data 

Generation Control Manipulator Conversion 


Figure 10. Overall Experiment Design. 

A MATLAB waveform generation software program was developed and used to 
produce discrete-time versions of selected continuous-time input signals to serve as input 
excitations. The input excitations were then transferred to an IBM compatible 386 PC 
where another software algorithm written in the ‘C’ programming language was used to 
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modulate various waveform parameters and to channel the input excitations to specified 
joints of the HMTB manipulator. While the HMTB manipulator arm responded to the 
input excitations, a data acquisition program written in the Ada programming language 
was used to extract joint data from the 1553 bus and to record the data to the HMTB 
Control Computer. This data includes the actual input/output time history from each joint. 
The recorded 1553 data file was then converted to an ASCII flat file format using an 
additional algorithm developed using the Matlab language. The ASCII flat file containing 
input/output time histories was then used for nonparametric and parametric analysis. 

Prior to the identification experiments, assumptions were made as to the model 
form and bandwidth of the open-loop dynamics for each joint. First, a PD control model 
with feedforward torque was assumed for each joint. This assumption was based on 
analytical models of the flight arm and partial documentation for the ground-based 
manipulator. Because of the size of the HMTB manipulator as well as its intended 
purpose, the bandwidth for each joint was assumed to correspond to astronaut response 
times (3 to 5 Hertz). Due to this assumption, all input excitations used in this 
identification were limited to 10 Hertz to satisfy the Nyquist requirement. At least a 
second order model was expected due to proportional (P) and derivative (D) components 
initially assumed for each PD control loop. The following sections will provide more 
detail on the experiment design segments. 

2, Input Excitations 

According to Soderstrom and Stoica [9], the input signal used in an identification 
experiment can significantly influence the resulting parameter estimates. Also, certain 
system identification methods require special types of inputs depending on the type of 
identification to perform. With this in mind, various input signals (excitations) were used 
to identify the dynamic parameters of the hydraulic manipulator. Many of the input 
excitations used in the system identification experiments were considered ‘normal’ test 
signals such as simple sinusoids, sum of sinusoids, pseudorandom binary sequences, and 
chirp input signals. A bipolar ramping pulse test signal was also used. This signal was 
used as a means of determining amplitude response characteristics of the joints in 
question. All test signals, however, were used as input excitations in determining single- 
input, single-output (SISO) as well as multi-input, multi-output (MIMO) black box 
models of the selected hydraulic joints. All input excitations were fed through the hand 
controller interface. The actual input/output data used for system identification, however, 
were extracted from the 1553 bus which recorded the measured and commanded angles at 
each joint. 

A. Description of Input Excitations 

This section will describe the input signals used to excite the HMTB manipulator 
for the purpose of system identification. The rationale for selecting the signals will also be 
discussed along with a general declaration of properties and characteristics for each input 
signal. 

As stated previously, several signals (input excitations) were used in determining 
the dynamic parameters of the dexterous orbital serving system manipulator arm. These 
input excitations include simple sinusoids, sum of sinusoids, pseudorandom binary 
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sequences, bipolar ramping pulses, and chirp signals. Each signal was used to excite three 
joints (shoulder yaw, shoulder pitch, and elbow pitch) of the HMTB arm. Modulation of 
signal parameters will be discussed later. 

Good identification experiments provide informative data by which different 
models can be discriminated within an intended model set. To provide this informative 
data, persistently exciting (pe) input signals must be selected [10], An input signal u(n) is 
said to be pe of order m if the spectral density <X>(w) is not equal to zero for at least m 
points in the interval - n < w < n for discrete-time systems. With the spectral density 
<D(w) defined as the discrete Fourier transform of the correlation function, that is, 

«>('*•) ■ T- Z R4r (2.1) 

r=-oo 

it has also been determined that u(n) is pe of order m if the limit of the autocorrelation 
function exists, that is. 


R u (r)= lim -J- ^u(t + r)u T (t) , 
and the autocorrelation matrix 


R«( n ) = 


R»(0) 

R u( -1) 


R u (l) 

R u (0) 


R u (n-1 ) 


R u (l-n) .. R u (0) 


( 2 . 2 ) 


(2.3) 


is nonsingular [10], The white noise signal eft), simulated in Figure 11, is persistently 
exciting of all orders [9], 

White Noise 
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Sinusoidal Input 

For the most simple input excitation, a sinusoid was selected to provide frequency 
response characteristics at a particular frequency and phase shift. Because of its 
simplicity, this signal was selected primarily for test purposes, that is, to determine if the 
joint in question was indeed responding to the input signal. Using a sinusoidal input 


u(t) = a sin( wt ), (2.4) 

the steady state output, assuming the system is linear, will become 

y(t) = b sin( wt ) (2.5) 

where 

b = a |G(hv)|, and (2.6) 

<t> = arg [G(iw)]. (2.7) 


The phase <f> for this signal will be negative, else the system is responding with no input. 
It should be noted that this input excitation is rather sensitive to disturbances (noise) and 
could be improved by repeating the sinusoid at a number of frequencies to obtain a 
graphical representation of the transfer function G(w) as a function of w. 

Sum of Sinusoids 

A sum of sinusoids provides a slight variation from the simple sinusoid by 
increasing the number of sinusoidal inputs with distinct frequency components which 
yields a greater bandwidth (BW) in the frequency domain. This type of input is used 
primarily in transfer function analysis. The discrete-time sum of two sinusoids input 
expressed mathematically as 

u(n) = ai sin( w/n ) + sin( wjh ) (2.8) 

where 

0 < W] < w 2 < 2 n (2.9) 

was used in the identification experiments, though the number of sinusoids need not be 
limited to two. As a general rule, the input u(n) will be pe of order 2m where m is the 
number of sinusoids in the sum [9], Therefore, to identify a fourth order system, only two 
sinusoids need be summed. The power spectral density for an infinite sequence of the sum 
of sinusoids is 


a 


a , 


O(w') = — [ S(w - w, ) + S(w + w, ) ] + — [<?(v»'-v»' 2 ) + <5(w' + w' 2 )]. (2.10) 
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Since there are exactly m=4 nonzero points in the interval the actual input signal 

u(n) is said to be pe of order at least 4. Figure 1 2 displays a sum of two sinusoids signal 
where aj = 1.0, a 2 = 1.0, wj = 0.02 ;r, and w 2 = 0.08 ;r. Figure 13 displays a plot of the 
estimated power spectral density of the sum of sinusoids signal. Plots of the actual 
excitation signals used for system identification of the HMTB joints are shown in the next 
chapter along with their power spectral densities. 

For notation purposes, the term ‘sum of sinusoids at 10 Hertz’ used in this thesis 
means that the highest frequency component f 2 of the given condition f 2 = 4 f\ for the sum 
of two sinusoids input will be equal to 10 Hz. This also implies that the lower frequency 
component// is equal to 2.5 Hz. For all sum of two sinusoids experiments performed on 
the HMTB joints, the arbitrarily selected condition f 2 = 4fi will hold. 




Pseudorandom Binary Sequences 

Due to their easy generation, the pseudorandom binary sequence (PRBS) has been 
a convenient input signal for many identification methods. The PRBS signal shifts 
between two levels in a certain pattern such that its first- and second-order characteristics, 
the mean value jj. and the correlation function R/ r ), are very similar to those of a white 
noise process e(t) provided that the number of samples used in the PRBS calculation is 
large. Interpretations vary as to the actual number of samples used, but is usually 
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experiment dependent. Therefore, this type of input signal applies itself well to determine 
correlation effects of various system parameters. Since the PRBS is band-limited and is 
periodic, it differs from a true white noise process. The PRBS is said to be pe of order 
equal to its period. Soderstrom and Stoica [9] indicate that in most cases the period of a 
PRBS is chosen to be of the same order as the number of samples in the experiment, or 
larger. Figure 14(a) depicts the PRBS signal. Its mathematical expression can be realized 
as 

u(t)= (C, + C 2 ) + (C,-C 2 )sign(R(r ) u(t-l) + w(t) ) (2.11) 

where 

Ci and C 2 are permissible binary levels, 

R(r ) is the covariance function, and 
w(t) is a simulated white noise process. 

Figure 14(b) displays the power spectral density of the PRBS input signal. 
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Chirp Signal 

Chirp is a technique invented by B. M. Oliver at Bell Labs in which signals are 
represented by a rapid up or down sweep in frequency [13], Chirp signals (sine sweeps) 
have been used for both radar and communication applications. This waveform was 
chosen as an input signal because of its selectable frequency range. Chirp signals have 
been known to produce regions with low power spectrum [14], For this reason, Franklin, 
Powell, and Workman [14] describe an expression for a chirp signal that does not have 
low power spectrum in the desired bandwidth. Their chirp waveform is expressed as 


r* = Ao + a*sin2;r f k k, 

(2.12) 

( k\ N -k\ 

a k = cw sat — — sat — — , and 
\pN) pN J 

(2.13) 

k 

fk ~ f start — {/stop ~f start ) > 

N 

(2.14) 

where 

N - number of points in the data window, 

Ao = constant reference offset adjustment, 
a max = maximum amplitude, 

p = fraction of window length for amplitude ramps, 
f start - starting frequency of chirp, and 
fuop = stopping frequency of chirp. 


The chirp expression used in the identification experiments of this thesis, however, can be 
characterized with the following equation: 

u(t) = cos( 2tt (J, +/A) + 27rf 0 t ) 

(2.15) 


where 

fi is the starting frequency, 
fh is the ending frequency, 
fo is the center frequency, and 
A =/*-//• 

A chirp waveform is shown in Figure 15(a) having values /; = 1 Hz,/ 0 = 4 Hz, and f h = 2 
Hz. The power spectral density of the chirp is shown in Figure 15(b). 

Bipolar Ramping Pulse 

The bipolar ramping pulse shown in Figure 16 has been included primarily to test 
the amplitude response characteristics of the manipulator joints. This input signal 
produces a series of periodic, alternating, ramping pulses. The user defined pulse width 
remains constant throughout the pulse sequence. In system identification literature, the 
impulse and step inputs have been used for transient analysis for several nonparametric 
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Figure 16. (a) Bipolar Ramping Pulse (b) PSD of Bipolar Ramping Pulse. 
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identification experiments. This sequence of alternating, ramping pulses will serve as a 
test signal for parametric model estimation. 


B . Generation of Input Excitations 

All input excitations were generated by a program written using MATLAB. The 
program allows a user to graphically display the generated waveform after selecting its 
appropriate parameters, such as its frequency and amplitude. The software then outputs a 
100-point waveform data file representing a discrete-time version of the continuous-time 
signal. Most data files contained one complete cycle of the waveform in order to allow 
the excitation control computer to accurately control the frequency of the selected 
discrete-time waveform. The length of the actual excitation is a periodic version of the 
100 samples. In this way, each excitation was allowed to reach steady-state conditions. 
Figure 17 shows the initial user interface for the waveform generation software. 



Figure 17. Waveform Generation Software User Interface. 

3, Input Excitations and Control 

As seen in the overall experiment design depicted in Figure 10, the generated 
waveform is controlled and channeled to various joints of the HMTB manipulator by the 
excitation control computer. This section will describe the hardware, procedures, and 
algorithms used to channel the generated excitation signals to specific joints of the HMTB 
manipulator for the purpose of conducting single-input, single-output (SISO) and multi- 
input, multi-output (MIMO) system identification tests. 


18 




A. Single Joint Excitation 

The excitation control computer, an EBM 386 PC, was used to perform all system 
identification tests. This PC contained the AT-MIO-16 National Instruments data 
acquisition board [15]. The AT-MIO-16 applies itself well to various multifunction 
analog, digital, and timing applications. In the experiment design, the AT-MIO-16 was 
used to convert the discrete-time input waveform to an equivalent frequency modulated 
analog signal. This was accomplished by forming a periodic version of each 100-sample 
waveform and then outputting the new frequency modulated signal through the digital-to- 
analog converter. 

Since the purpose of the identification experiment was to identify a dynamic model 
localized about an RPCM trajectory vector, the excitation control computer was 
interfaced to the HMTB manipulator at the exact location astronauts would be interfaced, 
that is, at the interface for the 2x3 DOF hand controller. Figure 18 shows the wiring 
implementation used to interface the AT-MIO-16 data acquisition board to the 2x3 DOF 
hand controller processor for SISO identification tests. 


AT-MIO-16 board 


2x3 DOF HC processor 



Figure 18. Single Axis Interface Configuration. 


Each single axis test consisted of moving the HMTB arm autonomously to the 
RPCM initial trajectory called the RPCM HOLSTER OUT APPROACH POINT. Next, 
the HMTB manipulator was placed in single joint mode allowing only the selected joint to 
accept input while all other joints were actively servoing. With the RPCM ORU (Figure 
9) loaded in the HMTB end effector, input control was transferred to the excitation 
control computer to perform the SISO tests in Table 1. All SISO tests were performed in 
position mode with direct input to each joint variable. 

Code written in the ‘C’ language was used to modify the control frequency and 
number of waveform iterations of the AT-MIO-16 board. For safety reasons, all AT- 
MIO-16 single axis input excitations were first tested on an oscilloscope. 
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Single Input/Single Output Tests 


Input Excitation 


single sinusoid 


sum of two sinusoids 


sum of two sinusoids 


sum of two sinusoids 


seudorandom binary sequence 


seudorandom binary sequence 


bipolar ramping pulse 


chirp signal 



Time/Frequency Characteristics 


freq=5 Hz 


freqi=1.25 Hz, ffeq2=5 Hz 


freqi=1.25 Hz, freq2=5 Hz 


freq j =2. 5 Hz, ffeq2=10 Hz 


order<100 


order<100 


freq= 1 Hz, pulse width=0. 1 sec 


freq start 5 Hz, freqend 10 Hz 


Table 1. Single-Input, Single-Output Tests. 

B, Multiple Joint Excitation with Bias Compensation 

Multi-input, multi-output (MIMO) tests, shown in Table 2, were used to identify 
the dynamical characteristics of the HMTB joints as well as to verify the models obtained 
using the SISO identification tests. The excitation control computer to HMTB 
manipulator wiring interface used in the SISO experiments was modified for the MIMO 


Multiple Input/Multiple Output Tests 

Test 

No. 

Input Excitation 

Maximum 

Amplitude 

Time/Frequency Characteristics 

1 

single sinusoid 

1 

freq =5 Hz 

2 

single sinusoid 

2 

freq=5 Hz 

3 

single sinusoid 

2 

freq=10 Hz 

4 

sum of two sinusoids 

1 

freq i=l. 25 Hz, freq2=5 Hz 

5 

sum of two sinusoids 

2 

ffeqi=1.25 Hz, freq2=5 Hz 

6 

sum of two sinusoids 

2 

freq i=2. 5 Hz, freq2=10 Hz 

7 

pseudorandom binary sequence 

1 

orderdOO 

8 

pseudorandom binary sequence 

2 

orderdOO 

9 

bipolar ramping pulse 

4 

freq=l Hz, pulse width=0.1 sec 

10 

chirp signal 

1 

freq s tart=0 Hz, freq en d=l Hz 

11 

chirp signal 

1 

freq s tart=0 Hz, freq en d=5 Hz 

12 

chirp signal 

1 

freq sl art = 5 Hz, freq en d=10 Hz 

13 

chirp signal 

1 

freq s iart=5 Hz, freqend=15 Hz 


Table 2. Multi-Input, Multi-Output Tests. 
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tests. This modification involved designing and implementing a bias compensator to offset 
the 2x3 DOF hand controller biases in the X, Y, and Z translational axes. Figure 19 shows 
the bias compensating wiring scheme used to interface the AT-MIO-16 data acquisition 
board to the 2x3 DOF hand controller processor for MEMO identification tests. 


AT-MIO-16 board 2x3 DOF HC processor 



Figure 19. Multiple Axes Interface with Bias Compensation. 

Preliminary procedures for MIMO testing involved moving the HMTB arm 
autonomously to the RPCM HOLSTER OUT APPROACH POINT and then transferring 
input control to the excitation control computer. The HMTB manipulator was then placed 
in Cartesian mode allowing only the translational inputs to accept values with respect to 
the end effector control frame. All rotational inputs (Euler angles) were held as constant 
as possible. With the RPCM ORU loaded in the HMTB end effector, the excitation 
control computer performed the MIMO tests in Table 2. All MIMO tests were performed 
in position mode. 

4. 1 553 Bus Data Acquisition 

As each identification test was performed on the HMTB manipulator, an Ada 
software program recorded various joint parameters. The 1553 data acquisition program 
recorded, measured, and commanded joint angles from sensors located at the manipulator 
actuators as well as force and moment data from the force-torque sensor located at the 
end effector. Data was recorded at 50 Hz, which is the fixed position loop transfer rate. 
This rate served to provide the Nyquist sampling frequency (25 Hz) for the input 
excitations used. In the identification experiments, the excitations were well below the 
Nyquist frequency. The Nyquist rate, however, was not as important as the bandwidth of 
the excitations. 

Another constraint imposed on the experiment design was the limited available 
memory storage for recording the joint data. Approximately twenty minutes of recording 
time was allotted on the control station computer. This equates to recording a total of 


21 



five experiment tests per run. After each set of five tests, the recorded files would have to 
be transferred to another computer to provide memory for another set of tests. 

5 . 1553 Bus Format to ASCII Conversion 

The final segment of the experiment design involved converting the 1553 recorded 
data file to an ASCII flat file format. This task was performed by a MATLAB conversion 
program. The conversion program extracted measured and commanded joint data from 
the 1553 formatted data file to be identified and saved this data in an ASCII flat file format 
to be used for nonparametric and parametric analysis. 


22 



m. Nonparamctric Model Estimation 


As a precursor to parametric identification, nonparametric methods are used first 
to extract characteristics of the unknown system which provides information in how to 
apply various parametric techniques. This chapter will show results of applying frequency, 
correlation, and spectral analysis techniques to the shoulder yaw, shoulder pitch, and 
elbow pitch joints of the HMTB manipulator. The results will help determine appropriate 
parametric model structures for the next chapter. 

1 . Procedure Description and Rationale 

Nonparametric model estimation involves determining a system’s characteristics 
from Bode plots and plots of input/output cross-correlation. Though sufficient, 
nonparametric methods give only moderately accurate models. For time domain 

nonparametric analysis, the impulse response and the step response are both useful in 
determining some basic control related characteristics of a system such as delay time, 
static gain and dominating time constants. Frequency domain techniques provide 
information such as the estimated transfer function, the bandwidth of a system, and a 
system’s phase characteristics. The techniques employed in this investigation include 
transfer function analysis, correlation analysis, and spectral analysis. 

Transfer function analysis was used to determine the frequency response of the yet 
to be identified system. This information helped to determine the frequency range of the 
input excitations to be used in the identification experiments. The frequency response 
approach was performed by applying a sum of sinusoidal inputs to the system and then 
recording the input/output time histories for each joint. Autocorrelation and cross- 
correlation functions were first computed from the data and then transformed to power 
spectral density and cross-power spectral density estimates, respectively. Spectral 
estimates were smoothed and averaged by using a Hamming window with the lag length 
approximately equal to a tenth of the number of data points. The estimated transfer 
function for each joint was computed as the ratio of the cross-power spectrum to the input 
power spectrum. Each joint’s transfer function estimate was represented in the form of 
Bode plots. 

Correlation analysis techniques were employed to provide information on the 
degree of linear dependence of a system’s parameters, that is, how well future values of 
the data can be predicted based on past observations. Correlation analysis is usually based 
on white noise or any input signal that is independent of the disturbances. A distinct 
advantage of correlation techniques is its insensitivity to additive noise on the output [9], 

Spectral analysis, a very versatile nonparametric technique, used various 
persistently exciting input signals to yield spectral estimates of the system. The spectral 
density or spectrum is a frequency domain function used to measure the frequency 
distribution of the mean square value of the data. Spectral estimates for each joint were 
computed using a Hamming window with the lag length approximately equal to a tenth of 
the number of data points. 
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2. Shoulder Yaw Joint 
A. Transfer Function Analysis 

Since astronauts will use the hand controller to operate the manipulator joints and 
the input signals for this identification were introduced through the same interface, the 
frequency range of the input signals were selected to coincide with astronaut response 
times (3 to 5 Hz) [16]. To perform transfer function analysis of the shoulder yaw joint, a 
sum of two sinusoids input whose frequency content ranged from 2.5 Hz to 10 Hz was 
introduced into the shoulder yaw position loop. The upper frequency (10 Hz) was 
selected because it was at least twice the average frequency response for astronauts. 
Deductively, if the identified model is valid for twice the intended bandwidth then it is 
reasonable to assume that the actual model will be well behaved within the intended 3 to 5 
Hz bandwidth. Figure 20 shows both sum of sinusoids input and output discrete time 
sequences recorded for the shoulder yaw joint during the identification experiment. As 
seen in Figure 20, considerable noise is present on the input signal. From the actual 
output sequence in this same figure, it is clear that the data is affected by disturbances. 
This is perhaps due to background interference being carried through the hardware 
interface. Figure 21 shows a magnified version of the shoulder yaw joint waveforms. This 
version of the input sequence shows the effect of the sample-and-hold and quantization 
functions being implemented by the HMTB control computer on the hand controller 
signal. 


2 
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Figure 20. Shoulder Yaw Sum of Sinusoids I/O at 10 Hertz. 
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Figure 2 1 . Shoulder Y aw Sum of Sinusoids Activity. 

If the dynamics of the shoulder yaw joint is assumed to be linear in a small 
localized region, then the output y( t ) can be seen as a weighted sequence of the form 

y(0 = S h(k)u( t-k) + v(t) (3.1) 

Jr=0 

where 

h(k) is the weighting sequence, and 
v(t) is the disturbance. 

The autocorrelation function may be estimated from the input data sequence as follows: 

Ru(t)= lim -J-£w(/ + r)u r (0. (3.2) 

Note that the cross-correlation function R^( r ) may be computed in the same manner. 
Taking the Fourier transform of 

Ryu(T ) = 2 h(k) R„(r - k ) (3.3) 

k = 0 

yields 

H(e iw )0„(w). (3.4) 
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The estimated transfer function is computed as 
H( e' iw ) = <&„(*) / 0». 


(3.5) 


The discrete-time transfer function estimate for the sum of sinusoids data sequence is 
shown in Figure 22. 
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Figure 22. Shoulder Yaw Transfer Function Estimate. 

Graphical interpretation of this transfer function yields a gain less than one for the 
entire bandwidth with a low frequency cutoff at approximately w= 6 rad/sec. The 
negative slope slightly above the break frequency indicates a second-order system until 
approximately 10 rad/sec. The rest of the graph shows additional resonances and 
disturbances of the system above 10 rad/sec. 

B. Correlation Analysis 

Correlation analysis techniques were applied to the shoulder yaw joint to provide 
information on the degree of linear dependence of the input and output of the joint. 
Correlation analysis is usually based on white noise or any input signal that is independent 
of the disturbances. For this reason, a pseudorandom binary sequence (PRBS) was used 
to excite the joint. PRBS signals simulate white noise statistical properties for the purpose 
of nonparametric identification. The one difference between PRBS and white noise is its 
periodicity. The mathematical PRBS expression has already been shown (2.1 1). Figure 
23 shows the entire input/output sequence of the PRBS input signal applied to the 
shoulder yaw joint during the identification experiment. Figure 24 shows a magnified 
version of the PRBS shoulder yaw joint activity. 



AMPLITUDE PLOT, input # 1 output # 1 
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Explanation of correlation analysis can be discussed by using the definition of the 
covariance function, that is, 

0,(0= E[{x(t). n x }{y(t+ r )- Hy}] (3.6) 

where 

E[ ] is the expectation operator, 

/i m is the expected value (mean) of sequence m. 

It can be shown that the covariance function and the correlation function are related 
through the following relationship, 

Cxy(r ) = Rxy( r ) - HxHy. (3.7) 

Since the PRBS has zero mean, the covariance and correlation functions are equivalent. 

Figure 25 shows three graphical representations of the output covariance (the 
autocorrelation of the output), the autocorrelation of the input, and the cross-correlation 
from the input to the output. The first graph in Figure 25 shows how the output signal is 
correlated with the transfer function. The autocorrelation of the input shows a signal that 
is white in nature but exhibits some periodicity as can be seen by the small graphical peaks 
which is expected for a PRBS. The autocorrelation graph of the input is typical since the 
autocorrelation function is an even function. The autocorrelation function evaluated at 
zero yields the mean square value of the input. The cross-correlation graph displays 
propagation characteristics of the joint such as the distance and/or the velocity of an input 
through the system. Cross-correlation also gives an estimate of the order of the system. 
The peaks of the cross-correlation graph indicate the contribution of each of several 
independent sources of excitation found in the output measurement. 

OUTPUT #1 
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Figure 23 . Shoulder Yaw PRBS Data. 
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Figure 24. Shoulder Yaw PRBS Activity. 
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C. Spectral Analysis 

Spectral analysis, the Fourier transform of the autocorrelation function, was used 
to measure the frequency distribution of the mean square value of the data. Two input 
excitations were used to determine the spectra of the joint. The bipolar ramping pulse 
(BRP) whose energy focused around 1 Hz was used while a sum of two sinusoids input 
was used with frequency components at 2.5 and 10 Hz. 

The bipolar ramping pulse signal was used to determine several spectral estimates 
of the shoulder yaw joint (Figure 26). Each spectral estimate was computed using a 
Hamming window. 


OUTPUT #1 

0.02r 
0 . 01 - 
0- 

- 0 . 01 - 
- 0 . 02 - 

.0 031 ' ‘ ‘ 1 1 ' 

0 5 10 15 20 25 30 


0. 


- 0 . 

0 5 10 15 20 25 30 


INPUT #1 


1 

. 1 

II 

■ 

1 1 1 

1 -* 



Figure 26. Shoulder Yaw Bipolar Ramping Pulse Data. 

Plots of the estimated disturbance spectrum, the output spectrum, the input power 
spectrum, and the cross-spectrum are shown in Figures 27 through 30, respectively. The 
plots indicate that disturbance phenomena are predominantly focused in the lower 
frequencies around one hertz. The output spectrum and the input power spectral density 
(PSD) both show higher amplitudes in the one to five hertz range. Various higher 
frequency lobes indicate that other modes of the system are being excited by the input 
signal. The estimated cross-spectrum reveals the same information. 

The small spectral amplitudes are attributed to the very small hand controller gains 
used in the HMTB control computer. That is, the HMTB computer system scaled the HC 
signals to a very small range before allowing the input signals to enter the position control 
loop. Since the experiment control computer, used to perform the system identification 
tests, was interfaced through the hand controller hardware, it was subject to the same 
input scaling resulting in small spectral amplitudes. 
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Estimated Disturbance Spectrum 



Figure 27. Shoulder Yaw Estimated Disturbance Spectrum. 



Figure 28. Shoulder Yaw Estimated Output Spectrum. 
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Power Spectral Density (PSD) 




Figure 30. Shoulder Yaw Estimated Cross-Spectrum. 

The sum of sinusoids input signal was used to determine several spectral estimates 
for frequencies less than 10 Hertz (62.8 rad/sec). Each spectral estimate was computed 
using a Hamming window. Figure 3 1 shows the entire I/O data record for the shoulder 
yaw joint. Figure 32 shows the sum of sinusoids activity. 
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Spectral estimates, Figures 33 through 36, reveal that the disturbances are at least 
two orders of magnitude lower for frequencies less than 10 Hertz (62.8 rad/sec). The 
estimated cross-spectrum reveals that other modes of the system are being excited. 

10 
5 
0 
-5 

0 5 10 15 20 25 

0. 

0 . 


0 5 10 15 20 25 

Figure 31. Shoulder Yaw Sum of Sinusoids Input. 
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Figure 32. Shoulder Yaw Sum of Sinusoids Activity. 
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Power Spectral Density (PSD) 



Figure 35. Shoulder Yaw Estimated Power Spectrum. 



frequency (rad/sec) 


Figure 36. Shoulder Yaw Estimated Cross-Spectrum. 
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3. Shoulder Pitch and Elbow Pitch Joints 

Identical nonparametric procedures performed for the shoulder yaw joint were also 
performed for the shoulder pitch and elbow pitch joints. Transfer function analysis for 
both joints yielded break frequencies at approximately seven radians per second. 
Correlation analysis as applied to both joints yielded two to five delay units from the input 
to the output of each joint implying possible system orders. Results from spectral analysis 
for each joint indicated that the disturbances were at least an order of magnitude lower 
than the output spectrums. Plots and graphs from the nonparametric procedures described 
above are displayed in Appendix A. 1 for both shoulder pitch and elbow pitch joints. 

4. Nonparametric Conclusions 

Transfer function analysis, correlation analysis, and spectral analysis techniques 
have been used to determine a crude nonparametric estimation of three HMTB 
manipulator joints (shoulder yaw, shoulder pitch, and elbow pitch). The nonparametric 
model estimation techniques used in this chapter suggest that parametric models should be 
selected to properly model the noise dynamics as well the system’s dynamics. 

Nonparametric analysis described each joint with minor to moderate process and 
measurement disturbances. The plots reveal greater measurement disturbances than 
process disturbances. Transfer function analysis of each joint indicates that models need 
to be constructed within a 1 to 5 Hz bandwidth. 

Errors in the nonparametric estimations may be attributed to several sources: 
random errors, bias errors, quantization errors in the experiment design, and choice of 
input signal. Random errors are caused by nonlinearities in the system. Bias errors are 
due to resolution errors in the spectral density estimates as well as unmeasured inputs that 
contribute to the output. Quantization errors are caused by the sample-and-hold function 
used in the HMTB control computer when accepting the hand controller input signal. 
Velocity limits in the control system also contributed to errors in the estimations. The 
experiment design introduced errors on the input measurements with improperly shielded 
wires in the experiment control computer interface. When these errors are introduced into 
an experiment design, the likelihood of obtaining accurate estimates decreases. 
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IV. Parametric Model Estimation: 

Transfer Function and State-Space System Identification 


The parametric identification procedures employed in this chapter use various 
black-box transfer function model structures to determine model estimates for the HMTB 
joints. The transfer function models, found within the MATLAB System Identification 
Toolbox by Lennart Ljung, use prediction error techniques to determine parameters for 
each black-box model. Residual analysis and cross-validation procedures will be primarily 
used to choose the best model estimate for each joint. 

1. Parametric Procedures 

Though sufficient, nonparametric methods as discussed in the previous chapter 
give only moderately accurate models from observed input/output data. To obtain more 
accurate model estimates, parametric identification techniques are used. The basic 
requirements for parametric identification are the observed input/output data, a set of 
candidate model structures, and a criterion to select the best model in the set [10], The 
system identification process as described by Ljung is shown in Figure 37, that is, after 
data has been collected from an experiment, a model structure is chosen, the criterion to 
identify a particular model in the structure is selected, the model is then calculated and 
validated. If the model is not satisfactory, another criterion is selected or another 
structure is chosen. Ljung’ s parametric identification process is quite iterative [10], 



Figure 37. System Identification Process. 

Model structures tested for the identification of the shoulder yaw, shoulder pitch, 
and elbow pitch joints include the autoregressive with extra input (ARX) model, the 
autoregressive moving average with extra input (ARMAX) model, the output error (OE) 
model, and the four-stage instrument variable (IV4) model forms. These model structures 
were selected to produce the best approximation for each joint’s dynamic characteristics. 
During the parameter estimation and analysis procedures, the ARX and IV4 structures did 


36 




not produce consistent results, therefore, only results from the ARMAX and OE model 
structures will be shown. These results are consistent with the results obtained from the 
nonparametric tests performed earlier. The nonparametric estimation yielded considerable 
information about the noise dynamics of each joint which coincides with the fact that both 
ARX and IV4 model structures do not sufficiently characterize the noise dynamics. 

For each HMTB joint, parametric techniques will be employed to determine the 
best model that fits several data sets. To perform this task, a transfer function that 
corresponds to the model will be obtained, residual analysis will be performed to 
determine the whiteness and independence of the model estimate’s equation errors, and 
pole-zero plots will be shown to determine if the model estimate is stable. The model will 
be compared to the I/O data to determine if the estimate produces a proper fit. Next, 
cross-validation will be shown to determine if the model estimate can fit other data sets. 
The state-space representation of the best model estimate will be obtained. And finally, 
the linear combination of state-space representations will be determined to produce a 
multivariable state-space estimate of all three 
joints. 

Figure 38 shows the operator interface developed specifically for this thesis to 
perform nonparametric estimation and parametric evaluation of the HMTB joints. The 
algorithms used in the evaluation code utilize functions from the MATLAB System 
Identification Toolbox (version 3). 




System Identification Toolbox Main Menu 



Figure 38. System Identification Operator Interface. 
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2. Parametric Black-Box Models 

Most n-th order systems can be described with a simple, linear difference equation 

y(t) + a, y(t- 1 )+...+ a„ a y(t-n a ) = 

b,u(t- 1 )+... + b„ b u(t-n b ) + £(t). (4.1) 

The disturbance term £( t ) serves as a direct error in the difference equation. This 
general model is generally referred to as an equation error model. The linear block-box 
models used in this section serve to estimate the general equation error model. The 
equation error dynamical model may also be described as 


y(t) ~G(q, O)u(t) + H(q, 0) £(t) 

where 

G(q ) is the system transfer function, 

H(q) is the disturbance transfer function, 

£(t) is the disturbance, 

0 is the parameter vector, and 

q is the delay operator. 

The ARMAX linear block-box model structure corresponds to setting 

G(<I> ~ <? "* and H(q) - ^ 

A(q) A(q ) 

where 


(42) 


(4.3) 


C(q) - 1 + Ci q' 1 + c 2 q* 2 + ... + Cnc q^ . 

The ARMAX structure gives considerable freedom in describing the properties of the 
disturbance term by estimating the error equation as a moving average of white noise. 
This structure describes the system that has a common factor in the denominators of the 
G(q) and H(q) polynomials. 

The output error (OE) model structure allows the transfer functions, G( q ) and 
H( q ), to be independently determined. That is, neither transfer function has a common 
polynomial description. The OE structure has the model form 

y( O = u( t) + £(t). (4.4) 


Both ARMAX and OE models structures are estimated using a prediction error 
method (PEM). The PEM is a modification of the least squares (LS) method. In the 
general LS method, the estimation procedure is performed by selecting the parameter 
vector 6 that minimizes the loss function, 
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(4.5) 


V(» ) - 

The PEM enhances the LS approach by forming the residual £ ( t ) as a difference 

A 

between the measured output y( t ) and a prefiltered output prediction y ( t\t-l; 6), that 
is, 

8(t) = y( t) - y ( t\t-l; 6) (4.6) 

where 


y ( t\t-l ; 6) = 6) (F V; 0)u(t) + { I - ITV; 0)}y(t). 

For a model estimate to correctly describe an unknown system, the residuals (equation 
errors) must be ideally white and independent of the input. 

3 , Identification of Shoulder Yaw Joint 
A. Preliminary Model Estimates 

Several ARX, ARMAX, IV4, and OE model structures were used to identify the 
dynamical characteristics of the shoulder yaw joint. Among these model structures, only 
the ARMAX and OE estimates exhibited a better fit among many data sets. Therefore, 
only ARMAX and OE estimates will be discussed. From the experiment tests, the 
shoulder yaw appeared to exhibit a more nonlinear response. This information quickly 
implies perhaps a higher order design to estimate this joint’s response. Throughout the 
identification process, models selected have been those that were the simplest to obtain 
while yet maintaining stability and the best approximation to many data sets. 

ARMAX 

After several iterations, an ARMAX model containing five poles, two zeros, and 
two delays on the input was found to sufficiently characterize the shoulder yaw joint. The 
transfer function is expressed as 

0.003804 z 3 + 0.003504 z 2 - 0.004285 z 

H( z ) (4 7) 

z 5 - 2.542 z 4 + 2.367 z 3 - 0.8523 z 2 - 0.08044 z + 0.1103 

Prediction error analysis of the ARMAX estimate yielded residuals that were white and 
with a high degree of independence. Figure 39 shows the residuals of the ARMAX model 
estimate using the PRBS input signal. The first plot in Figure 40 shows the whiteness of 
the model’s residuals (the autocorrelation of the residual) while the second (lower) plot 
shows the residual independence (the cross-correlation of the residual and the input) as a 
function of lag (delay). The dashed lines in each plot represent 99% confidence intervals. 
That is, if the curve in each plot goes significantly outside the confidence intervals, the 
model is not accepted as a good estimate. 
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Figure 39. Shoulder Yaw ARMAX Residuals. 



Figure 40. Shoulder Yaw ARMAX Residual Whiteness and Independence. 
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The pole-zero plot (Figure 41) displays the poles and zeros of the ARMAX model 
estimate. Since the poles of the discrete-time system are in the unit circle, the model is 
stable. The close pole and zero in the graph indicate a near pole-zero cancellation possibly 
indicating that the model order selected was too high. The other models tested without 
the close pole and zero produced estimates that did not sufficiently characterize the 
dynamics. 


OUTPUT#! INPUT#! 



- Figure 41 . Shoulder Yaw ARMAX Pole-Zero Plot. 

The ARMAX model estimate was then compared to the data set that produces the 
model. Figure 42 shows a comparison of the estimated model output to the measured 
output. Even though the model didn’t follow the PRBS data set very well, it showed the 
best flexibility in following many other data sets. Figure 43 shows the cross-validation of 
the ARMAX model estimate to the sum of sinusoids ten hertz input signal. Since cross- 
validation of a model estimate is a much harder task, the cross-validation of the model to 
various data sets weighed heavily in determining the best model estimate. 
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Output#! Fit: 0.01015 



Figure 42. Shoulder Yaw ARMAX Output Comparison. 



Figure 43. Shoulder Yaw ARMAX Cross-Validation. 
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OE 

An OE model containing two poles, one zero, and two delays on the input also 
produced good results for the shoulder yaw joint. The transfer function for this OE 
estimate is 


0.0464 z + 0.06144 

H( z ) = . (4.8) 

z 3 - 0.7501 z 2 - 0.1815 z 

Prediction error analysis of this OE estimate yielded residuals that were not very white. 
This might be attributed to the fact that the OE structure focuses more on the dynamics G 
and less on the noise properties H. Figure 44 shows the residuals of the OE model 
estimate using the PRBS input signal. The first plot in Figure 45 shows the whiteness of 
the model’s residuals (the autocorrelation of the residual) while the second (lower) plot 
shows the residual independence as a function of lag (delay). The cross-correlation of the 
OE residuals to the input signal shows slight independence for small positive and negative 
lags. Since the residual independence does not go significantly outside the confidence 
intervals, the OE model is accepted as a possible shoulder yaw estimate. 


Residuals of Current Estimate 



Figure 44. Shoulder Yaw OE Residuals. 

The pole-zero plot (Figure 46) displays the poles and zeros of the OE model 
estimate. Since the poles of the discrete-time system are inside the unit circle, the model is 
stable. Various tests using a higher number of OE poles yielded an unstable system. 
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Figure 45. Shoulder Yaw OE Residual Whiteness and Independence. 



Figure 46. Shoulder Yaw OE Pole-Zero Plot. 
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Output # 1 Fit: 0.004963 



Solid (-) : Model output, Dot (.) : Measured output 

Figure 47. Shoulder Yaw OE Output Comparison. 

. Q 3 Output # 1 Fit: 0.0006806 



Solid (-) : Model output, Dot (.) : Measured output 
Figure 48. Shoulder Yaw OE Cross-Validation. 

The OE model estimate was then compared to the data set that produces the 
model (Figure 47). The OE model estimate only fitted the mean of the PRBS data rather 
than the peaks of the sequence. This fit is perhaps better than the ARMAX comparison, 
however, the OE showed poor flexibility in following other data sets. Figure 48 shows the 
cross-validation of the OE model estimate to the sum of sinusoids ten hertz input signal. 
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The OE estimate exhibited the correct frequency response, as seen by the zero crossing in 
Figure 48, but failed to match the peaks of the sequence. This response is due to the 
magnitude of the OE estimate being less than unity at the frequency of the input data. 
This OE response was typical to other data files. 

B, Determination of Best Model Estimate 

The ARMAX structure with five poles, two zeros, and two delays on the input 
was selected as the best linear approximation for the shoulder yaw dynamics in a localized 
region. Since the bode plots of the ARMAX estimate (Figure 49) showed comparable 
results to the OE and ARX estimates shown, the ARMAX model exhibited some of the 
true characteristics of the unknown system. The ARMAX estimate was selected because 
of the whiteness and independence of its residuals. Also, the cross-validation of the model 
to several data sets described a more flexible estimate. For these reasons, the ARMAX 
model was selected as the best estimate of the shoulder yaw joint. 



Figure 49. Bode Plots of the ARMAX, ARX, and OE Estimates. 


C. SISO State-Space Estimate 

The SISO state-space estimate for the ARMAX shoulder yaw joint is: 

x(k + 1 ) = Aix(k) + Biu(k) + K, e(k) (4.9) 

y(k) = Cix(k) + Diu(k) + e(k) 


where 
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2.5420 1.0000 

-2.3667 0 
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4. Identification of Shoulder and Elbow Pitch Joints 

The same model structures (ARX, ARMAX, IV4, and OE) and procedures used to 
identify the dynamical characteristics of the shoulder yaw joint were also used to identify 
parametric dynamic models for the shoulder pitch and elbow pitch joints. Plots and graphs 
from the shoulder pitch and elbow pitch joints are shown in Appendix A.2. Summary of 
the best model estimate for each joint follows: 

Shoulder Pitch 

The OE model containing two poles, one zero, and one delay on the input was 
selected as the best linear approximation for the shoulder pitch dynamics in a localized 
region. The transfer function for this OE estimate is 

0.03459 z + 0.07584 

H( z ) = . (4.10) 

z 2 - 0.9776 z + 0.07922 

There are two primary reasons for selecting the OE estimate as the best dynamic 
approximation. First, the residuals were truly independent of the input which implies that 
the model is a very good approximation to the real joint dynamics. And second, the cross- 
validation of the model to several data sets described a more flexible model estimate. The 
SISO state-space estimate for the OE shoulder pitch joint is: 

x(k+l) = A 2 x(k) + B 2 u(k) + K 2 e(k) (4.11) 

y(k) = C 2 x(k)+D 2 u(k) + e(k) 

where 
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a 2 


0.9776 1.0000 

-0.0792 0 


B 2 


0.0346 

0.0758 


C 2 = [1.0000 0], 


D 2 = 0, 


K 2 = 


given X0 = 


Elbow Pitch 

An ARMAX model containing two poles, one zero, and two delays on the input 
was found to sufficiently characterize the dynamics of the elbow pitch joint. The transfer 
function for this model estimate is expressed as 

0.01237 z + 0.0426 

H( 2 ) = . (4.12) 

z 3 - 1.488 z 2 + 0.5435 z 

There are two primary reasons for selecting the ARMAX estimate as the best dynamic 
approximation. First, the residuals were whiter than the OE estimate indicating that the 
ARMAX appropriately modeled the noise characteristics of the joint. And second, the 
cross-validation of the model to several data sets described a more flexible model estimate. 
The SISO state-space estimate for the ARMAX elbow pitch joint is: 

x(k+J) = A 3 x(k) + B 3 u(k) + K 3 e(k) (4.13) 

y(k ) = C 3 x(k) + D 3 u(k) + e(k) 


where 



■ 1.4877 1.0000 0 


0 

a 3 = 

-0.5435 0 1.0000 

b 3 = 

0.0124 


0 0 0 


0.0426 

C 3 = [1.0000 0 0], d 3 

= o, 




' 1.4698 * 


'O' 

k 3 = 

-0.2739 

, given X0 = 

0 


0 


0 
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5. State-Space Multivariable Representation 

Each SISO state-space estimate previously determined was identified about an 
operating point. For the purpose of this identification, the operating point was chosen to 
coincide with the insertion point trajectory vector for the RPCM ORU exchange task. 
With each input signal producing only small deviations around the operating point, a local 
neighborhood was defined about the RPCM insertion point for which each identified SISO 
model is valid. This is the essence of linear approximation of a nonlinear model [18], 

The following three-joint, linear, multivariable state-space estimate was formed 
using each joint’s best dynamic SISO representation: 

x(k+ 1 ) = Ax(k) + Bu(k) + Ke(k) (4.14) 


y(k ) = Cx(k) + Du(k) + e(k) 


where 


A = 


A 

o 

0 


0 

A 

0 


0 

0 , 
A. 



0 

A 

o 


0 

0 , 
*3. 


c = 


C, 

0 

0 


0 

C 2 

0 


0 

0 

c. 


K = 


K, 

0 

0 


0 

K 2 

0 


0 

0 

K, 


D = 


D. 

0 

0 


0 0 
d 2 0 

0 D, 


The matrices A^, B„, C„, D„, and K n where n = 1, 2, or 3 refer to matrices previously 
determined in equations 4.9, 4.11, and 4.13. This multivariable state-space estimate will 
be used for comparison purposes. The actual RPCM data set along with several MIMO 
data sets will be used to cross-validate this multivariable estimate. 
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V. Parametric Model Estimation: 

Observer/Kalman Filter Identification 

Among alternate system identification procedures are the ones based on system 
realization theory. One such technique, used in the identification of space structures, is 
Observer/Kalman Filter Identification (OKID). This technique computes Markov 
parameters from pulse system response histories using an asymptotically stable observer to 
form a stable discrete state-space model. This chapter will briefly discuss the OKID 
technique provided in the System/Observer/Controller Identification Toolbox (SOCIT) by 
Jer-Nan Juang, Lucas G. Horta, and Minh Phan [19]. Residual analysis and cross- 
validation procedures will be used to identify the best state-space models for the HMTB 
joints. 


1 . OKID Background and Procedure 

When a pulse sequence is used as input into the discrete-time state-space dynamic 
equation 

x(k + 1 ) = A x(k) + Bu(lc) ( 5 . 1 ) 

y(k) = C x(k) + D u(k), 


the resulting series of equations can be formed into a pulse-response matrix F, that is, 

Y=[D CB CAB CA k, B ]. (5.2) 

The elements of the matrix F are known as the system Markov parameters. System 
realization involves determining the matrices A, B, C, and D from the system Markov 
parameters to satisfy the state and measurement equations (5.1) [8], Minimum realization 
theory, attributed to Ho and Kalman [12], determines a model with the smallest state- 
space dimension among all realizable systems. The procedure starts by forming the 
generalized Hankel matrix composed of Markov parameters in the following manner: 



' Y k 

n*. 

Y 1 

2 k+fi- 1 

H(k -1) = 

n*. 


y 

A k+p 



n + . • 

Y 

1 k+a+fl-2 


(5.3) 


If the number of rows a and the number of columns /? are greater than the order of the 

system, then the Hankel matrix is of rank n. This realization algorithm extracts linear 
state-space matrix components from noise-free data. 

For noisy input/output sequences, the Eigensystem Realization Algorithm (ERA) 
produces better results [8], By deleting specific rows and columns of the Hankel matrix, 
ERA forms a block data matrix composed of strongly measured Markov sequence 
components. A minimum realization may be obtained by factorization of the block data 
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matrix using singular value decomposition. The order of the identified system is 
determined by selecting the number of significant singular values. 

Observer/Kalman Filter Identification is determined by inserting an observer into 
the discrete-time state-space dynamic equation (5.1) to form the discrete-time state-space 
observer model 

x(k + J ) = Ax(k)+Bv(k) (5.4) 

y(k) = Cx(k) + Du(k), 


where 

A = A + GC, 


B = [ B+GD -G J, and 


v(k) = 


u{k) 

y(k)_ ‘ 


When a pulse sequence is used as input into the observer model (5.4), the following 
observer Markov parameters may be computed: 

Y — [ D CB CAB CA k, B ]. (5.5) 

The OKID technique then computes the system Markov parameters from the observer 
Markov parameters for minimum realization of a state-space model estimate. It is obvious 
from (5.5) that the matrices A, B, C, and D are embedded in the observer Markov 
parameters. Since the observer gain G may be arbitrarily chosen, the OKID routine 

creates a deadbeat observer by simply placing all the eigenvalues of A at the origin 
producing an asymptotically stable observer. Setting G to be the deadbeat observer gain 
allows for a minimum number of Markov parameters to describe the input/output 
relationship of a system [8], 

Juang [8] describes the relationship between the state-space observer model and 
the Kalman filter equation 

X(k + J) = AX(k) + Bu(k) + K£ r (k) (5.6) 


y (k) = CX(k) + Du(k), 

where 

/v 

X(k) is the estimated state, 

A 

y (k) is the estimated measurement, 
K is the Kalman filter gain, and 
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£ r (k) is the residual defined as the difference y(k) -y (k) . 

The observer gain is said to be the steady-state Kalman filter gain 

K = -G (5.7) 

in theory if the residuals are identically zero, S r ( k ) = 0, and the data length is 

sufficiently long to produce negligible truncation error. Theoretical background of the 
OKID technique is found in the text Applied System Identijicgtign by Jer-Nan Juang [8], 

For each joint, system and observer parameters will first be determined. Next, the 
associated prediction errors will be computed. The Hankel matrix will be shown for 
proper order selection. After selecting the system order from the Hankel singular values, 
the state-space estimate will be realized. This realization will also yield the Kalman filter 
gain which for the purposes of this investigation will be approximated to the observer gain 
since each estimate will be selected to minimize the residuals and the arm will be operated 
in a localized region to minimize system nonlinearities. Each model estimate will be 
compared to the data set that produced the model as well as cross-validated to another 
data set. A three-joint multivariable state-space model will be determined from the three 
SISO state-space estimates. 

2, Identification of Shoulder Yaw Joint 
A. Determine Markov Parameter Set 

An upper bound for the OKID model order must first be specified to compute an 
estimate. Using a fifth order system as the upper bound, five independent observer 
parameters were initially computed to identify the shoulder yaw state-space model using 
the PRBS input/output data. The system and observer Markov parameters for the 
shoulder yaw joint are shown in Figure 50. As seen in Figure 50, the rate of decay for the 
observer parameters is much larger than that for the system Markov parameters. This 
demonstrates the advantage of the deadbeat observer in minimizing the number of pulse 
response samples used to realize the state-space equation [8]. The relatively high number 
of independent Markov parameters shown in the plot indicates that the shoulder yaw joint 
exhibits relatively significant nonlinear characteristics. 
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Pulse Response Estimate of System 



Pulse Response Estimate of Observer 



Figure 50. Shoulder Yaw System and Observer Markov Parameters. 

B. OKID State-Space Estimate 

The normalized prediction errors associated with the independent observer 
parameters and the input data file are shown in Figure 5 1 . The lower plot in Figure 5 1 
shows the variance of each observer parameter with the measured data. A smoothing 
error is also plotted next to each variance. Using the Hankel matrix of singular values 
(Figure 52), a system order of two was selected for minimum realization of the shoulder 
yaw joint, that is, the second-order model obtained described 99.8691% of the test data. 
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Figure 5 1 Shoulder Yaw Output Prediction Errors. 



The following discrete-time state-space observer model has been realized for the 
shoulder yaw joint: 

x(k+ 1) = (A, + G,Ci)x(k) + [ Bj+GjDj - G,Jv(k ) (5.8) 

y(k) = Cix(k) + D,u(k), 


54 






where 


A, = 


.99487 .02896 

-.15355 55941 


2.7065 

12.366 


Ci = [004706 -.001125], D, = 0.000714, 


G, = 


-276.41 

-123.20 


and v( k ) 


C. Analysis 

The bode plots of the discrete-time state-space model estimate is shown in Figure 
53. The frequency range of this estimate seems to be valid for extremely low frequencies 
less than 1 Hz. 

Response due to input= 1 output= 1 
50 r • — ■ " ' r>l — ' — I — 1 — ' ' ' —■'l 



10 10 
Frequency (Hz) 



10 10 10 10 1CT 

Frequency (Hz) 

Figure 53. Shoulder Yaw Bode Plots of State-Space Estimate. 


To validate the state-space model estimate, the model output was compared to the 
data set that produced the model. Figure 54 shows the predicted state-space output 
compared to the measured output. As seen in the graph, the second-order model estimate 
predicts the measured output extremely well. 
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Output Number = 1 



Time (sec) [ (.) Measured Output, (-) Estimated Output ] 
Figure 54. Shoulder Yaw Predicted Versus Measured Output. 


Cross-validating the shoulder yaw model OKID estimate with the five to ten Hertz 
linearly swept chirp signal (Figure 55) produced good results. Residuals for the cross- 
validation were very low. This implies that the deadbeat observer can possibly be 
considered the Kalman Filter gain according to 5.7. 


Output Number = 1 



Time (sec) [ (.) Measured Output, (-) Estimated Output ] 


Figure 55. Shoulder Yaw Chirp Cross-Validation. 
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It should be noted that numerous iterations were performed to obtain this state- 
space model estimate. In the iterations where the data was not detrended, the best 
estimate for minimum realization was a fourth-order system. This higher-order system 
produced extremely poor residuals and therefore was rejected as the shoulder yaw 
estimate. The only advantage of the fourth-order system was the minimum number of 
independent observer Markov parameters needed to realize the state-space dynamic 
equation. 

3 . Identification of Shoulder and Elbow Pitch Joints 

The same OKID procedures used to identify the dynamic characteristics of the 
shoulder yaw joint were also used to identify parametric dynamic models for the shoulder 
pitch and elbow pitch joints. Plots and graphs from the shoulder pitch and elbow pitch 
joints are shown in Appendix A.3. Summary of both model estimates follows: 


Shoulder Pitch 

A system order of two was selected for minimum realization of the shoulder pitch 
joint. More specifically, the second-order model obtained described 96. 1075 % of the test 
data. The following discrete-time state-space observer model has been realized for the 
shoulder pitch joint: 


x(k + 1 ) — (A 2 + G 2 C 2 ) x( k) + [ B 2 +G 2 D 2 - G 2 ] v(k) 
y(k) = C 2 x(k) + D 2 u(k), 


where 


A 2 


.90891 .10801 

-.35206 .33925 


B 2 


6.881 
11.762 ’ 


C 2 = [.014742 -.0065875], 


G2 - 


-65.283 

22.597 


and 


D 2 = 
v(k) 


0.010405, 

'«(*)■ 
y(k\ ■ 


(5.9) 


Elbow Pitch 

A system order of two was selected for minimum realization of the elbow pitch 
joint. More specifically, the second-order model obtained described 98.8832 % of the test 
data. The following discrete-time state-space observer model has been realized for the 
elbow pitch joint: 

x( k + 1 ) = (A3 + G3C3) x( k ) + [ B3+G3D3 - G 3 J v( k ) (510) 

y(k) = C 3 x(k) + D 3 u(k), 

where 
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A 3 = 


.94738 .077488 

-.27646 .517280 ’ 


Bs 


4.7407 
11.0670 ’ 


Cs = [.013258 -.0059914], D 3 = 0.0032941, 


G 3 


-77.049 
11.586 ' 


and 


v(k) = 


u(k) 
y{k)_ ' 


Though higher order models may have produced better fits in reducing the residuals, it 
was more advantageous to minimize the system order thus reducing the complexity of 
both state-space dynamic estimates. 

4, State-Space Multivariable Representation 

Each SISO state-space estimate previously determined was identified about the 
insertion point for the RPCM ORU exchange task. With each input signal producing only 
small deviations around the operating point, a local neighborhood was defined about the 
RPCM insertion point for which each identified SISO model is valid. The following three- 
joint, linear, multivariable state-space estimate was formed using the dynamic SISO state- 
space representation for each joint: 

x(k + 1 ) = Ax(k) + Bu(k) + Ke(k) ( 5 . 11 ) 

y(k) = Cx(k) + Du(k) + e(k) 

where 



X 

0 

0 “ 


X 

0 

o' 





A = 

0 

A 2 

0 

B = 

0 

B 2 

0 

J 





0 

0 

Aj _ 


0 

0 

B 3 _ 






"c, 

0 

0 " 


X 

0 

0 

-1 


X 

0 

0 ' 

C = 

0 

C 2 

0 

D = 

0 

d 2 

0 


K = 

0 

k 2 

0 


_ 0 

0 

C 3. 


0 

0 

Ds 



0 

0 

K 3. 


The matrices A*, B n , C n , D„, and K„ where n = 1, 2, or 3 refer to matrices previously 
determined in equations 5.8, 5.9, and 5.10. This multivariable state-space estimate will be 
used for comparison purposes. The actual RPCM data set along with several MIMO data 
sets will be used to cross-validate this multivariable estimate. 
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VI. Comparison of Identified Models 


This chapter will compare the multivariable state-space model estimate obtained 
using prediction error techniques within Lennart Ljung’s System Identification Toolbox 
and the multivariable state-space model estimate obtained using the Observer/Kalman 
Filter identification (OKID) technique provided in the System/Observer/Controller 
Identification Toolbox (SOCIT) by Jer-Nan Juang, Lucas G. Horta, and Minh Phan. 
Frequency plots and pole/zero maps for each estimate will be shown. Both multivariable 
estimates will be compared to data sets obtained from an RPCM experiment and a MIMO 
experiment using the chirp signal. Fit comparisons and residual analysis will be performed 
for each state-space estimate. 

1 , Identified Model Forms 

The identification techniques investigated and described in the previous chapters 
represent parametric models of the form 

y(t) = G(q)u(t) + v(t) (6.1) 


where 

G(q) is the open-loop transfer function, and 
v(t) represents the disturbances. 

This linear equation attempts to describe the open-loop dynamic characteristics of each of 
the three HMTB joints with additive disturbances. It should be noted that this dynamic 
equation is, in essence, an open-loop representation of the actual closed-loop dynamical 
implementation for each joint. This implies that the actual closed-loop dynamics will be 
embedded within the open-loop description of each joint. For high noise-to-signal ratios in 
the I/O time histories of each joint, the disturbances may be represented as 

v(t) = H(q)e(t) (6.2) 

where 

H(q) is the disturbance dynamics, and 
e( t) represents white noise. 

In the System Identification Toolbox , prediction error techniques were used to 
determine ARMAX and OE black-box model structures for each joint. The ARMAX 
structure forms the joint dynamics according to the block diagram shown in Figure 56. 
The OE structure differs from the ARMAX structure by allowing the disturbances (noise) 
to go unfiltered as shown in Figure 57. The parametric black-box model developed for 
each joint was then converted to a SISO state-space estimate. 
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e(t) 



Figure 56. ARMAX Structure Block Diagram. 


e(t) 



Figure 57. OE Structure Block Diagram. 


A linear combination of the SISO state-space estimates were used to develop the 
multivariable state-space model. Linearization of the HMTB arm was performed by 
allowing only small perturbations in each joint’s input signal to produce a localized region. 
The multivariable estimate obtained using prediction error techniques described the 
dynamics in the flexible innovations discrete-time state-space model form 

x(k+ 1) = Ax(k) + Bu(k) +Ke(k) (6.3) 

y(k) = Cx(k)+Du(k) + e(k). 

The Observer/Kalman Filter Identification (OKID) technique identifies the dynamic 
characteristics of each joint using a discrete-time observer model of the form 

X (k + 1 ) = A X(k)+Bu(k) - G£ (k) (6.4) 

y(k) = C X(k) + Du(k) + £ (k). 

where 

X (k) is the estimate of state x (k), and 
£(k)=y(k) - y(k). 

It should be noted that the observer G can only be equated to the steady-state 
Kalman Filter gain K if and only if the residuals £ (k) are white, zero mean and Gaussian 
[8], In the following comparisons, the identified observer is not the Kalman Filter gain. 
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The computed observer G simply minimizes the residuals due to nonlinearities in each 
joint, non-whiteness of the noise processes, and correlation effects between the residual 
and the input signal. 

Both multivariable model estimates will be compared using equivalent state-space 
representations, that is K - -G where K in this case is simply a residual filter. For 
notation purposes, the multivariable model identified using prediction error techniques will 
be called the SysID model and the model identified using the Observer/Kalman Filter 
Identification (OKID) technique will be called the OKID model. 

2, Transfer Function Analysis 

The transfer function G(q) from equation 6.1 may expressed in terms of the state- 
space matrices A, B, C, and D as 

G(q) = C(ql - A)' 1 B + D. (6.5) 

Bode plots for the SysID MIMO model shown in Figures 58 - 60 describe the frequency 
and phase response characteristics of the shoulder yaw, shoulder pitch, and elbow pitch 
joints, respectively. 



Figure 58. SysID Shoulder Yaw Bode Plots. 
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Bode plots for the OKED MIMO model shown in Figures 61-63 describe the 
frequency and phase response characteristics of the shoulder yaw, shoulder pitch, and 
elbow pitch joints, respectively. As seen in the bode plots (Figures 58 - 63), the OKED 
model and the SysID model produced comparable results indicating that both techniques 
captured true frequency characteristics of each HMTB joint. 


62 

















Figure 63 . OKID Elbow Pitch Bode Plots. 


3. Pole/Zero Map 

Figure 64 shows the eigenvalues and transmission zeros of the SysID multivariable 
model. The eigenvalues (roots of the characteristic equation) represent the poles of the 
discrete-time MEMO estimate. Since the magnitude of the eigenvalues are less than unity, 
the discrete-time model is considered stable. The eigenvalue-transmission zero 
combination near the origin indicates that the SysID model may be unnecessarily complex 
to describe the HMTB joints. 
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Pole/Zero Map (X-Eigenvalues.O-Transmission Zeros) 



Real Axis 


Figure 64. SysID Multivariable Model Pole/Zero Map. 

The pole/zero map showing the eigenvalues and transmission zeros of the OKID 
multivariable model is shown in Figure 65. Since the magnitude of the eigenvalues in 
Figure 65 are less than unity, this model is considered stable. The OKID MIMO model 
differs from the SysID MIMO model in that the eigenvalues of the OKID model are all 
real, all transmission zeros are complex conjugates, and a minimum number of eigenvalues 
are used to describe the system. 



Real Axis 

Figure 65. OKID Multivariable Model Pole/Zero Map. 
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4, RPCM Fit Comparison 

Both multivariable state-space model estimates were compared to shoulder yaw, 
shoulder pitch, and elbow pitch data obtained from an RPCM experiment (joint angles 
shown in Figure 66) performed in the hydraulic manipulator test bed at NASA LaRC. 
Both estimates will be evaluated on the how well each model fits the data set as well as the 
whiteness and independence of each model’s residuals. Ideally, all residuals should be 
white and independent of the input for the model to perfectly identify the dynamic 
characteristics of the joints. 



Figure 66. Multivariable RPCM Experiment Data. 

SvsID RPCM Model Fit 

Figures 67 - 69 show outputs of the multivariable SysID model compared to actual 
outputs from the RPCM experiment. As seen in the Figures 67 and 69, the SysED model 
matches both the shoulder yaw and elbow pitch responses quite well. The shoulder pitch 
model output shown in Figure 68 is adequate but doesn’t quite match the peaks of the 
RPCM data set. This implies that the SysED model may not sufficiently characterize the 
true dynamic characteristics of the shoulder pitch joint. 
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Output Number = 1 



Figure 67. SysID RPCM Shoulder Yaw Comparison. 


Output Number = 2 



Time (sec) [ (.) Measured Output, (-) Estimated Output ] 
Figure 68. SysID RPCM Shoulder Pitch Comparison. 
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Figure 69. SysID RPCM Elbow Pitch Comparison. 

The residuals associated with the RPCM data and the SysID multivariable model 
estimate are shown in Figure 70. For the SysID model to correctly describe the dynamics 
of each joint, residuals must be ideally white and independent of the input. Figures 71 and 
72 show the whiteness of the SysID residuals associated with the RPCM data set. As 
seen in the figures, the shoulder yaw and elbow pitch residuals are fairly white. Residuals 
for the shoulder pitch joint, however, were not white. This is partially due to the use of 
the OE structure which focuses more on the dynamics than the noise properties of the 
shoulder pitch joint. 
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Correlation function of residuals. Output # 1 




Figure 71. SysID Shoulder Yaw and Pitch Residual Whiteness. 


Correlation function of residuals. Output # 3 



Cross corr. function between inpTjt 1 and residuals from output 1 



Figure 72. SysED Elbow Pitch Whiteness and Shoulder Yaw Independence. 

To determine residual independence for the multivariable SysID model estimate, 
the cross-correlation of the residual and the input for each joint was determined. The 
lower plot in Figure 72 shows the independence of the shoulder yaw residuals. For small 
positive lags, the shoulder yaw residuals are correlated. For lags greater than fifteen, the 
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SysID estimate correctly models the dynamics of the shoulder yaw joint for the RPCM 
data set. 


OKID RPCM Model Fit 

Figures 73 - 75 show outputs of the multivariable OKID model compared to actual 
outputs from the RPCM experiment. As seen in the comparison plots, the OKID model 
matches the shoulder yaw, shoulder pitch, and elbow pitch outputs very well. 


Output Number = 1 



Time (sec) [ (.) Measured Output, (-) Estimated Output ] 
Figure 73. OKID RPCM Shoulder Yaw Comparison. 
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Output Number = 2 
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Figure 74. OKID RPCM Shoulder Pitch Comparison. 
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Figure 75. OKID RPCM Elbow Pitch Comparison. 

The residuals associated with the RPCM data and the OKID multivariable model 
estimate are shown in Figure 76. The computed residuals for the shoulder yaw and 
shoulder pitch joints are fairly white as seen in Figures 77. Residuals for the elbow pitch 
joint (Figure 78, top plot) are not as white as the other joints. Overall, the OKID 
produces good residual whiteness for the RPCM experiment. 
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Figure 76. OKID RPCM M3MO Residuals. 


Correlation function of residuals. Output # 1 
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Correlation function of residuals. Output # 2 
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Figure 77. OKED Shoulder Yaw and Pitch Residual Whiteness. 
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Correlation function of residuals. Output # 3 



Cross corr. function between input 1 and residuals from output 1 



Figure 78. OK TP Elbow Pitch Whiteness and Shoulder Yaw Independence. 

The lower plot in Figure 78 shows the independence of the shoulder yaw residuals. 
For positive lags between five and fifteen, the shoulder yaw residuals were moderately 
correlated. This can possible be attributed to the significant nonlinearities previously 
found in the shoulder yaw joint. 

5, MIMO Chirp Fit Comparison 

Besides the RPCM experiment data, both multivariable state-space estimates were 
compared to various MIMO data sets to determine the constraints of each identified 
model. Results of comparing both model estimates to various MIMO data sets will be 
provided in the next section. This section will describe comparisons of the model 
estimates to the five to ten hertz linearly swept MIMO chirp data set (Figure 79). Output 
comparisons of each joint will be shown for both model estimates as well as the magnitude 
of each model’s residuals to the MIMO chirp data set. 
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Figure 79. Multivariable Chirp Experiment Data. 

SvsID MIMQ Chirp Fit 

Figures 80 - 82 show outputs of the multivariable SysID model compared to actual 
outputs from the MIMO chirp experiment. The SysID model effectively matched both the 
shoulder yaw and elbow pitch responses. The shoulder pitch model output (Figure 81) 
deviated slightly from the shoulder pitch chirp output. This characteristic was also found 
when the shoulder pitch model output was compared to the RPCM data set. This implies 
that the multivariable SysID model, though adequate, may not sufficiently characterize the 
true dynamics of the shoulder pitch joint. 
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Output Number = 1 



Time (sec) [ (.) Measured Output, (-) Estimated Output ] 

Figure 80. SysED Chirp Shoulder Yaw Comparison. 
Output Number = 2 



Time (sec) [ (.) Measured Output, (-) Estimated Output ] 
Figure 8 1 . SysID Chirp Shoulder Pitch Comparison. 
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Time (sec) [ (.) Measured Output, (-) Estimated Output ] 
Figure 82. SysID Chirp Elbow Pitch Comparison. 


The residuals associated with the MIMO chirp data set and the SysID multivariable 
model estimate are shown in Figure 83. The largest residual in the plot is attributed to the 
model errors of the shoulder pitch joint. Residual analysis for the MIMO chirp data set 
(not shown), were comparable to the results obtained using the RPCM experiment data. 
In other words, shoulder yaw and elbow pitch residuals are fairly white. Residuals for the 
shoulder pitch joint, however, were not white. 



Figure 83 . SysID MIMO Chirp Residuals. 
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Results of residual independence for the multivariable SysID model estimate 
yielded slight dependence for the shoulder yaw joint at small lags indicating that the SysID 
model may not adequately characterize shoulder yaw characteristics in the ten hertz range. 
A small amount of coupling was also found to exist between the shoulder yaw residuals 
and shoulder pitch input at this frequency range. 


OKID MIMO Chirp Fit 

Figures 84 - 86 show outputs of the multivariable OKID model compared to actual 
outputs from the MIMO chirp data set. The OKID model output matched the shoulder 
yaw, shoulder pitch, and elbow pitch outputs very well. 

Output Number = 1 
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Output Number = 2 



Time (sec) [ (.) Measured Output, (-) Estimated Output ] 


Figure 85. OKID Chirp Shoulder Pitch Comparison. 


Output Number = 3 



Time (sec) [ (.) Measured Output, (-) Estimated Output ] 


Figure 86. OKID Chirp Elbow Pitch Comparison. 


The residuals associated with the MIMO chirp data set and the OKID 
multivariable model estimate are shown in Figure 87. The OKID model estimate produced 
good residual whiteness for the shoulder pitch and the elbow pitch joints (not shown). 
However, shoulder yaw residuals showed only marginal whiteness. 
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Figure 87. OKID MIMO Chirp Residuals. 

Results of residual independence for the multivariable OKED model estimate 
produced significant correlation for the shoulder yaw joint at small lags indicating that the 
OKID model may not adequately characterize shoulder yaw dynamics in the ten hertz 
range. A small amount of coupling was found to exist between the shoulder yaw 
residuals and shoulder pitch input and between the shoulder yaw residuals and the elbow 
pitch input in the ten hertz frequency range. 

6, Comparison Results 

This section will compare both the SysID and OKID multivariable model estimates 
for several performance categories. For each category, the strengths and weaknesses of 
each model will be evaluated as well as the technique used to identify the model estimate. 
Performance categories will include frequency bandwidth, model stability, flexibility, 
parsimony, robustness, RPCM experiment fit, and various MIMO data fits. 

Frequency Bandwidth 

The SysID and OKID models both produced comparable Bode plots indicating 
that both techniques captured true frequency content of each HMTB joint. Both model 
estimates were able to follow data sets containing frequency content in the range of 3 to 5 
Hz though the identified cutoff frequency for each joint was found to be approximately 1 
Hz. Initially, this 1 Hz cutoff result was considered questionably low. Later, however, 
this low frequency cutoff was verified by a report conducted by the STX Corporation for 
NASA [16], The STX report found that Martin Marietta’s 1 Hz bandwidth was 
significantly lower than bandwidths recommended by prior research and by other research 
engineers [16], The uncertainty in proper bandwidth selection lies in the fact that an 
optimum bandwidth for space telerobot applications is unknown. 
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When compared to data sets containing a greater than 5 Hz frequency content, the 
OKID estimate slightly outperformed the SysID estimate for each joint. This may have 
been attributed to inadequate modeling of the SysID estimate. 

Stability 

Since the magnitude of the eigenvalues for both multivariable model estimates 
were less than unity, both models were found to be stable. 

Flexibility 

In terms of flexibility, the prediction error techniques used to determine the SysID 
model estimate provided numerous approaches to model both the dynamics of the system 
as well as the noise properties of the HMTB joints. The OKID technique, though more 
accurate, was less flexible. 

Parsimony 

The Observer/Kalman filter identification (OKID) technique was more 
parsimonious in its attempt to describe the dynamic characteristics of the HMTB joints. 
Information extracting using this technique produced a minimum realization in allowing 
each joint to be described by a second-order system. In contrast, the SysID model used a 
fifth-order model to describe the shoulder yaw dynamics, a second-order model to 
describe the shoulder pitch dynamics, and a third-order model to describe the elbow pitch 
dynamics. 

Robustness of estimates 

The SysID model proved to be more robust when compared to the OKID model 
estimate. Comparable SysID models were obtained even in the presence of bias errors and 
outliners in the data set. For instance, data detrending and/or filtering reduced a particular 
SysID model from a third-order system to a second-order system. In contrast, the OKID 
estimate was more sensitive to bias errors and outliners in the data set. One OKID model 
required a fourth-order system to describe the dynamics of a joint. However, a second- 
order model was sufficient when the data set was detrended. 

RPCM Experiment Fit 

The OKID multivariable model estimate provided a much better fit when 
compared to data obtained from the RPCM experiment. The OKID estimate obtained a 
more accurate fit by effectively minimizing the residuals for each joint model. The 
identified model showed minimal coupling between the joints in the localized region. The 
SysID model estimate consistently showed high residuals for the shoulder pitch SISO 
estimate. This can be attributed to the OE model structure used to characterize the joint. 

MIMO Data Fits 

The OKID estimate produced a better fit when compared to a majority of MIMO 
data sets. For low frequency data sets, the SysID estimate modeled the dynamics of the 
shoulder yaw joint much better than the OKID estimate. This would imply perhaps a 
hybrid model between the OKED estimate (shoulder pitch and elbow pitch) and the SysID 
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estimate (shoulder yaw) to provide an overall better multivariable model estimate for low 
frequencies. 

Overall, the Observer/Kalman filter identification (OKDD) technique produced a better 
multivariable estimate when compared to the SysID multivariable model estimate. Menu- 
driven software was developed and used to evaluate the OKDD model and to compare 
both multivariable models. 
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VIL Conclusions 


Linear, dynamic models for three joints of the hydraulic manipulator test bed 
(HMTB) at the NASA Langley Research Center have been identified using nonparametric 
and parametric system identification techniques. Nonparametric techniques yielded an 
approximate 1 Hz bandwidth for each joint using transfer function analysis, an expected 
order for each joint using correlation analysis, and the degree of process and measurement 
disturbances for each joint using spectral analysis. Two different parametric identification 
techniques were used and compared in developing dynamic models of the joints. The first 
parametric technique, used primarily for control system identification, employed a 
prediction-error method to produce a stable model estimate. The bandwidth for this 
estimate proved adequate when compared to several data sets. An advantage of this 
technique is its flexibility of use. The user has several options, alternatives, and methods 
from which to conduct an identification investigation. When compared to the RPCM 
experiment data, this technique yielded adequate results. The second parametric 
technique, used primarily in modal system identification, employed a minimum realization 
algorithm to produce a stable model estimate using only second-order systems to describe 
the characteristics of each joint. This technique was extremely simple to use while yet 
providing an adequate bandwidth for the identified models among many data sets. The 
models identified using this technique produced an accurate fit to both the RPCM 
experiment data and various MIMO data sets. 

Matlab menu-driven system identification software scripts were developed for this 
thesis. One program, using functions from the MATLAB System Identification Toolbox, 
was used to perform nonparametric and parametric analysis of the manipulator data The 
other program identified models using the Observer/Kalman Filter Identification (OKID) 
technique, provided in the System/Observer/Controller Identification Toolbox (SOCIT). 
The latter program used several toolboxes to perform MIMO comparisons for the 
identified multivariable models. Both programs were found to be extremely useful 
especially in minimizing the time required to perform nonparametric and parametric 
analysis of the data. The programs, modular in design, are easily expandable. Though 
written to use data from the manipulator joints, the programs may be easily adapted to 
incorporate data from any dynamic system for the purpose of identification. 

1 . Suggestions for Future Work 

With the identified models in this investigation valid only in a localized region with 
a specific loading, additional identification tests should be performed. Recursive 
identification methods, in particular, should be performed with several different loading 
configurations for the manipulator. Also, a more complete model should be identified by 
determining the dynamic parameters of the shoulder roll, wrist pitch, wrist yaw, and wrist 
roll joints (which were not identified in this thesis). 

2. Model Reference Control 

A model reference control system is proposed for the previously identified 
multivariable dynamic HMTB model. In this controls approach, the behavior of the 
ground-based model would closely resemble the behavior of the on-orbit flight model for 
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each joint. This capability would allow astronauts to perform crucial mission training 
tasks with a ground-based hydraulic manipulator that would be kinematically and 
dynamically equivalent to the flight manipulator. Figure 88 shows a block diagram of the 
model reference control system proposed for the DOSS manipulator. A distinct advantage 
of this control system is its ability to perform acceptably in the presence of nonlinearities, 
uncertainties, and variations in the identified system parameters [18], This service would 
de-emphasize errors developed in the dynamic parameters during the model identification 
process. 

Using the previously identified linear, dynamic HMTB model (equation 5.1), a 
linear model reference system can be described by the state equation 

x d (k+J) = Ax d (k) + Bv(k) (7.1) 


where 

Xd(k) is the state vector of the on-orbit dynamic model, 
v(k) is the input vector, 

A is the model reference state matrix, and 
B is the model reference state matrix. 

The model reference control system will have a stable equilibrium state if the magnitude of 
the discrete-time eigenvalues of A are less than unity. The primary task of the controller in 
Figure 88 will be to adjust the actuating signal of the identified HMTB model for the 
purpose of driving the state error between both models to zero [18], 



Figure 88. Model Reference Control System for DOSS. 
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Appendix A. Plots and Graphs 


This appendix contains plots and graphs for the shoulder pitch and elbow pitch 
joints. Explanation of the graphs are discussed in the main text. The graphs and plots are 
divided into three subdivisions as follows: 

A. 1 Nonparametric Plots 

A. 2 Classical Parametric Plots 

A. 3 OKID Parametric Plots. 
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A. 1 Nonparametric Plots 
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Figure A1 . Shoulder Pitch I/O Data. 
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Figure A2. Shoulder Pitch Transfer Function Estimate 
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Figure A3 . Shoulder Pitch PRBS Data. 
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Figure A4. Shoulder Pitch Correlation Plots. 
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Figure A5. Shoulder Pitch Bipolar Ramping Pulse Data. 
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Figure A6. Shoulder Pitch BRP Activity. 
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Estimated Output Spectrum 



Figure A1 3 . Shoulder Pitch Estimated Output Spectrum. 



Figure A 14. Shoulder Pitch Estimated Input Power Spectrum. 


92 







93 






AMPLITUDE PLOT, input # 1 output it 1 



frequency (rad/sec) 

Figure A17. Elbow Pitch Transfer Function Estimate. 
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Figure A1 8. Elbow Pitch PRBS Data. 
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Figure A20. Elbow Pitch Bipolar Ramping Pulse Data. 
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Figure A21 . Elbow Pitch Estimated Disturbance Spectrum. 
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Figure A22. Elbow Pitch Estimated Output Spectrum. 
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Figure A23. Elbow Pitch Estimated Input Power Spectrum. 



Figure A24. Elbow Pitch Estimated Cross-Spectrum. 
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Classical Parametric Plots 



Figure A30. Shoulder Pitch ARMAX Residuals. 



Figure A3 1 . Shoulder Pitch ARMAX Residual Whiteness and Independence. 
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Figure A32. Shoulder Pitch ARMAX Pole-Zero Plot. 
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Figure A3 3. Shoulder Pitch ARMAX Output Comparison. 
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Correlation function of residuals. Output # 1 
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Figure A38. Shoulder Pitch OE Output Comparison. 


Output # 1 Fit: 0.002869 



Figure A39. Shoulder Pitch OE Cross-Validation. 
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Figure A40. Bode Plots of the ARMAX, ARX, and OE Estimates. 



Figure A41 . Elbow Pitch ARMAX Residuals. 
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Correlation function of residuals. Output # 1 
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Figure A42. Elbow Pitch ARMAX Residual Whiteness and Independence. 
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Figure A43. Elbow Pitch ARMAX Pole-Zero Plot. 
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Figure A44. Elbow Pitch ARMAX Output Comparison. 
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Figure A45. Elbow Pitch ARMAX Cross-Validation. 
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Output # 1 Fit: 0.004003 
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Figure A46. Elbow Pitch OE Residuals. 
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Figure A47. Elbow Pitch OE Residual Whiteness and Independence. 
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Figure A48. Elbow Pitch OE Pole-Zero Plot. 


Output # 1 Fit: 0.003899 
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Figure A49. Elbow Pitch OE Output Comparison. 
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A. 3 OKID Parametric Plots 


Pulse Response Estimate of System 




Figure A52. Shoulder Pitch System and Observer Markov Parameters. 
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Figure A53. Shoulder Pitch Output Prediction Errors. 
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Figure A56. Shoulder Pitch Predicted Versus Measured Output. 
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Figure A57. Shoulder Pitch Chirp Cross-Validation. 
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Pulse Response Estimate of System 




Figure A58. Elbow Pitch System and Observer Markov Parameters. 
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Hankel Matrix Singular Values 
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Figure A60. Elbow Pitch Hankel Matrix. 
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Figure A61 . Elbow Pitch Bode Plots of State-Space Estimate. 
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Figure A62. Elbow Pitch Predicted Versus Measured Output. 
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Figure A63. Elbow Pitch Chirp Cross-Validation. 
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